
# Load packages
  library(tidyverse)
  library(foreign)
  library(lmtest)
  library(sandwich)
  library(stargazer)

# Calculate robust confidence intervals
  se_robust <- function(x)
    coeftest(x, vcov = vcovHC(x, type = "HC0"))[, "Std. Error"]
  p_robust <- function(x)
    coeftest(x, vcov = vcovHC(x, type = "HC0"))[, "Pr(>|t|)"]

  #### Dutch Parliamentary Election Study 2006
  
  # Load data
  dt <- read.spss("ObsStudy_00_data_Netherlands_2006.sav", to.data.frame=TRUE)
  
  # Party rating
  colnames(dt)[seq(237, 242, 1)] <- c("rat_PvdA_factor", "rat_VVD_factor", "rat_D66_factor", "rat_CDA_factor", "rat_SP_factor", "rat_GroenLinks_factor")
  
  # Generate numeric party rating variables
  dt$rat_PvdA <- ifelse(dt$rat_PvdA_factor == "0 unsympathetic", 0, NA)
  dt$rat_PvdA <- ifelse(dt$rat_PvdA_factor == "one", 1, dt$rat_PvdA)
  dt$rat_PvdA <- ifelse(dt$rat_PvdA_factor == "two", 2, dt$rat_PvdA)
  dt$rat_PvdA <- ifelse(dt$rat_PvdA_factor == "three", 3, dt$rat_PvdA)
  dt$rat_PvdA <- ifelse(dt$rat_PvdA_factor == "four", 4, dt$rat_PvdA)
  dt$rat_PvdA <- ifelse(dt$rat_PvdA_factor == "five", 5, dt$rat_PvdA)
  dt$rat_PvdA <- ifelse(dt$rat_PvdA_factor == "six", 6, dt$rat_PvdA)
  dt$rat_PvdA <- ifelse(dt$rat_PvdA_factor == "seven", 7, dt$rat_PvdA)
  dt$rat_PvdA <- ifelse(dt$rat_PvdA_factor == "eight", 8, dt$rat_PvdA)
  dt$rat_PvdA <- ifelse(dt$rat_PvdA_factor == "nine", 9, dt$rat_PvdA)
  dt$rat_PvdA <- ifelse(dt$rat_PvdA_factor == "10 very sympathetic", 10, dt$rat_PvdA)
  
  dt$rat_CDA <- ifelse(dt$rat_CDA_factor == "0 unsympathetic", 0, NA)
  dt$rat_CDA <- ifelse(dt$rat_CDA_factor == "one", 1, dt$rat_CDA)
  dt$rat_CDA <- ifelse(dt$rat_CDA_factor == "two", 2, dt$rat_CDA)
  dt$rat_CDA <- ifelse(dt$rat_CDA_factor == "three", 3, dt$rat_CDA)
  dt$rat_CDA <- ifelse(dt$rat_CDA_factor == "four", 4, dt$rat_CDA)
  dt$rat_CDA <- ifelse(dt$rat_CDA_factor == "five", 5, dt$rat_CDA)
  dt$rat_CDA <- ifelse(dt$rat_CDA_factor == "six", 6, dt$rat_CDA)
  dt$rat_CDA <- ifelse(dt$rat_CDA_factor == "seven", 7, dt$rat_CDA)
  dt$rat_CDA <- ifelse(dt$rat_CDA_factor == "eight", 8, dt$rat_CDA)
  dt$rat_CDA <- ifelse(dt$rat_CDA_factor == "nine", 9, dt$rat_CDA)
  dt$rat_CDA <- ifelse(dt$rat_CDA_factor == "10 very sympathetic", 10, dt$rat_CDA)
  
  dt$rat_VVD <- ifelse(dt$rat_VVD_factor == "0 unsympathetic", 0, NA)
  dt$rat_VVD <- ifelse(dt$rat_VVD_factor == "one", 1, dt$rat_VVD)
  dt$rat_VVD <- ifelse(dt$rat_VVD_factor == "two", 2, dt$rat_VVD)
  dt$rat_VVD <- ifelse(dt$rat_VVD_factor == "three", 3, dt$rat_VVD)
  dt$rat_VVD <- ifelse(dt$rat_VVD_factor == "four", 4, dt$rat_VVD)
  dt$rat_VVD <- ifelse(dt$rat_VVD_factor == "five", 5, dt$rat_VVD)
  dt$rat_VVD <- ifelse(dt$rat_VVD_factor == "six", 6, dt$rat_VVD)
  dt$rat_VVD <- ifelse(dt$rat_VVD_factor == "seven", 7, dt$rat_VVD)
  dt$rat_VVD <- ifelse(dt$rat_VVD_factor == "eight", 8, dt$rat_VVD)
  dt$rat_VVD <- ifelse(dt$rat_VVD_factor == "nine", 9, dt$rat_VVD)
  dt$rat_VVD <- ifelse(dt$rat_VVD_factor == "10 very sympathetic", 10, dt$rat_VVD)
  
  dt$rat_GroenLinks <- ifelse(dt$rat_GroenLinks_factor == "0 unsympathetic", 0, NA)
  dt$rat_GroenLinks <- ifelse(dt$rat_GroenLinks_factor == "one", 1, dt$rat_GroenLinks)
  dt$rat_GroenLinks <- ifelse(dt$rat_GroenLinks_factor == "two", 2, dt$rat_GroenLinks)
  dt$rat_GroenLinks <- ifelse(dt$rat_GroenLinks_factor == "three", 3, dt$rat_GroenLinks)
  dt$rat_GroenLinks <- ifelse(dt$rat_GroenLinks_factor == "four", 4, dt$rat_GroenLinks)
  dt$rat_GroenLinks <- ifelse(dt$rat_GroenLinks_factor == "five", 5, dt$rat_GroenLinks)
  dt$rat_GroenLinks <- ifelse(dt$rat_GroenLinks_factor == "six", 6, dt$rat_GroenLinks)
  dt$rat_GroenLinks <- ifelse(dt$rat_GroenLinks_factor == "seven", 7, dt$rat_GroenLinks)
  dt$rat_GroenLinks <- ifelse(dt$rat_GroenLinks_factor == "eight", 8, dt$rat_GroenLinks)
  dt$rat_GroenLinks <- ifelse(dt$rat_GroenLinks_factor == "nine", 9, dt$rat_GroenLinks)
  dt$rat_GroenLinks <- ifelse(dt$rat_GroenLinks_factor == "10 very sympathetic", 10, dt$rat_GroenLinks)
  
  dt$rat_SP <- ifelse(dt$rat_SP_factor == "0 unsympathetic", 0, NA)
  dt$rat_SP <- ifelse(dt$rat_SP_factor == "one", 1, dt$rat_SP)
  dt$rat_SP <- ifelse(dt$rat_SP_factor == "two", 2, dt$rat_SP)
  dt$rat_SP <- ifelse(dt$rat_SP_factor == "three", 3, dt$rat_SP)
  dt$rat_SP <- ifelse(dt$rat_SP_factor == "four", 4, dt$rat_SP)
  dt$rat_SP <- ifelse(dt$rat_SP_factor == "five", 5, dt$rat_SP)
  dt$rat_SP <- ifelse(dt$rat_SP_factor == "six", 6, dt$rat_SP)
  dt$rat_SP <- ifelse(dt$rat_SP_factor == "seven", 7, dt$rat_SP)
  dt$rat_SP <- ifelse(dt$rat_SP_factor == "eight", 8, dt$rat_SP)
  dt$rat_SP <- ifelse(dt$rat_SP_factor == "nine", 9, dt$rat_SP)
  dt$rat_SP <- ifelse(dt$rat_SP_factor == "10 very sympathetic", 10, dt$rat_SP)
  
  
  # Coalition evaluations/preferred coalition
  colnames(dt)[seq(229, 232, 1)] <- c("rat_coal_CDA_VVD_factor", "rat_coal_CDA_PvdA_factor", "rat_coal_PvdA_SP_GroenLinks_factor", "rat_coal_PvdA_VVD_factor")
  
  # Generate numeric coalition evaluation variable
  dt$rat_coal_CDA_VVD <- ifelse(dt$rat_coal_CDA_VVD_factor == "1 very undesirable", 1, NA)
  dt$rat_coal_CDA_VVD <- ifelse(dt$rat_coal_CDA_VVD_factor == "two", 2, dt$rat_coal_CDA_VVD)
  dt$rat_coal_CDA_VVD <- ifelse(dt$rat_coal_CDA_VVD_factor == "three", 3, dt$rat_coal_CDA_VVD)
  dt$rat_coal_CDA_VVD <- ifelse(dt$rat_coal_CDA_VVD_factor == "four", 4, dt$rat_coal_CDA_VVD)
  dt$rat_coal_CDA_VVD <- ifelse(dt$rat_coal_CDA_VVD_factor == "five", 5, dt$rat_coal_CDA_VVD)
  dt$rat_coal_CDA_VVD <- ifelse(dt$rat_coal_CDA_VVD_factor == "six", 6, dt$rat_coal_CDA_VVD)
  dt$rat_coal_CDA_VVD <- ifelse(dt$rat_coal_CDA_VVD_factor == "7 very desirable", 7, dt$rat_coal_CDA_VVD)
  
  dt$rat_coal_CDA_PvdA <- ifelse(dt$rat_coal_CDA_PvdA_factor == "1 very undesirable", 1, NA)
  dt$rat_coal_CDA_PvdA <- ifelse(dt$rat_coal_CDA_PvdA_factor == "two", 2, dt$rat_coal_CDA_PvdA)
  dt$rat_coal_CDA_PvdA <- ifelse(dt$rat_coal_CDA_PvdA_factor == "three", 3, dt$rat_coal_CDA_PvdA)
  dt$rat_coal_CDA_PvdA <- ifelse(dt$rat_coal_CDA_PvdA_factor == "four", 4, dt$rat_coal_CDA_PvdA)
  dt$rat_coal_CDA_PvdA <- ifelse(dt$rat_coal_CDA_PvdA_factor == "five", 5, dt$rat_coal_CDA_PvdA)
  dt$rat_coal_CDA_PvdA <- ifelse(dt$rat_coal_CDA_PvdA_factor == "six", 6, dt$rat_coal_CDA_PvdA)
  dt$rat_coal_CDA_PvdA <- ifelse(dt$rat_coal_CDA_PvdA_factor == "7 very desirable", 7, dt$rat_coal_CDA_PvdA)
  
  dt$rat_coal_PvdA_SP_GroenLinks <- ifelse(dt$rat_coal_PvdA_SP_GroenLinks_factor == "1 very undesirable", 1, NA)
  dt$rat_coal_PvdA_SP_GroenLinks <- ifelse(dt$rat_coal_PvdA_SP_GroenLinks_factor == "two", 2, dt$rat_coal_PvdA_SP_GroenLinks)
  dt$rat_coal_PvdA_SP_GroenLinks <- ifelse(dt$rat_coal_PvdA_SP_GroenLinks_factor == "three", 3, dt$rat_coal_PvdA_SP_GroenLinks)
  dt$rat_coal_PvdA_SP_GroenLinks <- ifelse(dt$rat_coal_PvdA_SP_GroenLinks_factor == "four", 4, dt$rat_coal_PvdA_SP_GroenLinks)
  dt$rat_coal_PvdA_SP_GroenLinks <- ifelse(dt$rat_coal_PvdA_SP_GroenLinks_factor == "five", 5, dt$rat_coal_PvdA_SP_GroenLinks)
  dt$rat_coal_PvdA_SP_GroenLinks <- ifelse(dt$rat_coal_PvdA_SP_GroenLinks_factor == "six", 6, dt$rat_coal_PvdA_SP_GroenLinks)
  dt$rat_coal_PvdA_SP_GroenLinks <- ifelse(dt$rat_coal_PvdA_SP_GroenLinks_factor == "7 very desirable", 7, dt$rat_coal_PvdA_SP_GroenLinks)
  
  dt$rat_coal_PvdA_VVD <- ifelse(dt$rat_coal_PvdA_VVD_factor == "1 very undesirable", 1, NA)
  dt$rat_coal_PvdA_VVD <- ifelse(dt$rat_coal_PvdA_VVD_factor == "two", 2, dt$rat_coal_PvdA_VVD)
  dt$rat_coal_PvdA_VVD <- ifelse(dt$rat_coal_PvdA_VVD_factor == "three", 3, dt$rat_coal_PvdA_VVD)
  dt$rat_coal_PvdA_VVD <- ifelse(dt$rat_coal_PvdA_VVD_factor == "four", 4, dt$rat_coal_PvdA_VVD)
  dt$rat_coal_PvdA_VVD <- ifelse(dt$rat_coal_PvdA_VVD_factor == "five", 5, dt$rat_coal_PvdA_VVD)
  dt$rat_coal_PvdA_VVD <- ifelse(dt$rat_coal_PvdA_VVD_factor == "six", 6, dt$rat_coal_PvdA_VVD)
  dt$rat_coal_PvdA_VVD <- ifelse(dt$rat_coal_PvdA_VVD_factor == "7 very desirable", 7, dt$rat_coal_PvdA_VVD)
  
  
  # Coalition likelihood
  colnames(dt)[seq(233, 236, 1)] <- c("coallik_CDA_VVD_factor", "coallik_CDA_PvdA_factor", "coallik_PvdA_SP_GroenLinks_factor", "coallik_PvdA_VVD_factor")
  
  # Generate numeric coalition likelihood variable
  dt$coallik_CDA_VVD <- ifelse(dt$coallik_CDA_VVD_factor == "1 very unlikely", 1, NA)
  dt$coallik_CDA_VVD <- ifelse(dt$coallik_CDA_VVD_factor == "two", 2, dt$coallik_CDA_VVD)
  dt$coallik_CDA_VVD <- ifelse(dt$coallik_CDA_VVD_factor == "three", 3, dt$coallik_CDA_VVD)
  dt$coallik_CDA_VVD <- ifelse(dt$coallik_CDA_VVD_factor == "four", 4, dt$coallik_CDA_VVD)
  dt$coallik_CDA_VVD <- ifelse(dt$coallik_CDA_VVD_factor == "five", 5, dt$coallik_CDA_VVD)
  dt$coallik_CDA_VVD <- ifelse(dt$coallik_CDA_VVD_factor == "six", 6, dt$coallik_CDA_VVD)
  dt$coallik_CDA_VVD <- ifelse(dt$coallik_CDA_VVD_factor == "seven", 7, dt$coallik_CDA_VVD)
  dt$coallik_CDA_VVD <- ifelse(dt$coallik_CDA_VVD_factor == "eight", 8, dt$coallik_CDA_VVD)
  dt$coallik_CDA_VVD <- ifelse(dt$coallik_CDA_VVD_factor == "nine", 9, dt$coallik_CDA_VVD)
  dt$coallik_CDA_VVD <- ifelse(dt$coallik_CDA_VVD_factor == "10 very likely", 10, dt$coallik_CDA_VVD)
  
  dt$coallik_CDA_PvdA <- ifelse(dt$coallik_CDA_PvdA_factor == "1 very unlikely", 1, NA)
  dt$coallik_CDA_PvdA <- ifelse(dt$coallik_CDA_PvdA_factor == "two", 2, dt$coallik_CDA_PvdA)
  dt$coallik_CDA_PvdA <- ifelse(dt$coallik_CDA_PvdA_factor == "three", 3, dt$coallik_CDA_PvdA)
  dt$coallik_CDA_PvdA <- ifelse(dt$coallik_CDA_PvdA_factor == "four", 4, dt$coallik_CDA_PvdA)
  dt$coallik_CDA_PvdA <- ifelse(dt$coallik_CDA_PvdA_factor == "five", 5, dt$coallik_CDA_PvdA)
  dt$coallik_CDA_PvdA <- ifelse(dt$coallik_CDA_PvdA_factor == "six", 6, dt$coallik_CDA_PvdA)
  dt$coallik_CDA_PvdA <- ifelse(dt$coallik_CDA_PvdA_factor == "seven", 7, dt$coallik_CDA_PvdA)
  dt$coallik_CDA_PvdA <- ifelse(dt$coallik_CDA_PvdA_factor == "eight", 8, dt$coallik_CDA_PvdA)
  dt$coallik_CDA_PvdA <- ifelse(dt$coallik_CDA_PvdA_factor == "nine", 9, dt$coallik_CDA_PvdA)
  dt$coallik_CDA_PvdA <- ifelse(dt$coallik_CDA_PvdA_factor == "10 very likely", 10, dt$coallik_CDA_PvdA)
  
  dt$coallik_PvdA_SP_GroenLinks <- ifelse(dt$coallik_PvdA_SP_GroenLinks_factor == "1 very unlikely", 1, NA)
  dt$coallik_PvdA_SP_GroenLinks <- ifelse(dt$coallik_PvdA_SP_GroenLinks_factor == "two", 2, dt$coallik_PvdA_SP_GroenLinks)
  dt$coallik_PvdA_SP_GroenLinks <- ifelse(dt$coallik_PvdA_SP_GroenLinks_factor == "three", 3, dt$coallik_PvdA_SP_GroenLinks)
  dt$coallik_PvdA_SP_GroenLinks <- ifelse(dt$coallik_PvdA_SP_GroenLinks_factor == "four", 4, dt$coallik_PvdA_SP_GroenLinks)
  dt$coallik_PvdA_SP_GroenLinks <- ifelse(dt$coallik_PvdA_SP_GroenLinks_factor == "five", 5, dt$coallik_PvdA_SP_GroenLinks)
  dt$coallik_PvdA_SP_GroenLinks <- ifelse(dt$coallik_PvdA_SP_GroenLinks_factor == "six", 6, dt$coallik_PvdA_SP_GroenLinks)
  dt$coallik_PvdA_SP_GroenLinks <- ifelse(dt$coallik_PvdA_SP_GroenLinks_factor == "seven", 7, dt$coallik_PvdA_SP_GroenLinks)
  dt$coallik_PvdA_SP_GroenLinks <- ifelse(dt$coallik_PvdA_SP_GroenLinks_factor == "eight", 8, dt$coallik_PvdA_SP_GroenLinks)
  dt$coallik_PvdA_SP_GroenLinks <- ifelse(dt$coallik_PvdA_SP_GroenLinks_factor == "nine", 9, dt$coallik_PvdA_SP_GroenLinks)
  dt$coallik_PvdA_SP_GroenLinks <- ifelse(dt$coallik_PvdA_SP_GroenLinks_factor == "10 very likely", 10, dt$coallik_PvdA_SP_GroenLinks)
  
  dt$coallik_PvdA_VVD <- ifelse(dt$coallik_PvdA_VVD_factor == "1 very unlikely", 1, NA)
  dt$coallik_PvdA_VVD <- ifelse(dt$coallik_PvdA_VVD_factor == "two", 2, dt$coallik_PvdA_VVD)
  dt$coallik_PvdA_VVD <- ifelse(dt$coallik_PvdA_VVD_factor == "three", 3, dt$coallik_PvdA_VVD)
  dt$coallik_PvdA_VVD <- ifelse(dt$coallik_PvdA_VVD_factor == "four", 4, dt$coallik_PvdA_VVD)
  dt$coallik_PvdA_VVD <- ifelse(dt$coallik_PvdA_VVD_factor == "five", 5, dt$coallik_PvdA_VVD)
  dt$coallik_PvdA_VVD <- ifelse(dt$coallik_PvdA_VVD_factor == "six", 6, dt$coallik_PvdA_VVD)
  dt$coallik_PvdA_VVD <- ifelse(dt$coallik_PvdA_VVD_factor == "seven", 7, dt$coallik_PvdA_VVD)
  dt$coallik_PvdA_VVD <- ifelse(dt$coallik_PvdA_VVD_factor == "eight", 8, dt$coallik_PvdA_VVD)
  dt$coallik_PvdA_VVD <- ifelse(dt$coallik_PvdA_VVD_factor == "nine", 9, dt$coallik_PvdA_VVD)
  dt$coallik_PvdA_VVD <- ifelse(dt$coallik_PvdA_VVD_factor == "10 very likely", 10, dt$coallik_PvdA_VVD)
  
  # Dependent variable: binary vote choice
  # v081 (V11_2) Vote intention [current] parliamentary elections
  dt$V081[dt$V081 %in% c("_duplicated_16", "DK yet", "blanc", "invalid")] <- NA
  
  # Generate party choice variable
  dt$vote <- dt$V081
  dt$vc_PvdA <- ifelse(dt$vote == "PvdA", 1, 0)
  dt$vc_CDA <- ifelse(dt$vote == "CDA", 1, 0)
  dt$vc_VVD <- ifelse(dt$vote == "VVD", 1, 0)
  
  # Coalition Evaluation
  Z <- dt %>% select("rat_coal_CDA_VVD", "rat_coal_CDA_PvdA", "rat_coal_PvdA_SP_GroenLinks", "rat_coal_PvdA_VVD") %>%
    mutate_all(as.numeric) %>%
    mutate_all(funs(scales::rescale(.,to = c(0, 1))))
  
  rat_coal_PvdA <- Z %>% select(rat_coal_CDA_PvdA, rat_coal_PvdA_SP_GroenLinks, rat_coal_PvdA_VVD) %>% as.matrix()
  rat_coal_CDA <- Z %>% select(rat_coal_CDA_VVD, rat_coal_CDA_PvdA) %>% as.matrix()
  rat_coal_VVD <- Z %>% select(rat_coal_CDA_VVD, rat_coal_PvdA_VVD) %>% as.matrix()
  rat_coal_GroenLinks <- Z %>% select(rat_coal_PvdA_SP_GroenLinks) %>% as.matrix()
  rat_coal_SP <- Z %>% select(rat_coal_PvdA_SP_GroenLinks) %>% as.matrix()
  
  # Coalition likelihood
  gamma <- dt %>% select("coallik_CDA_VVD", "coallik_CDA_PvdA", "coallik_PvdA_SP_GroenLinks", "coallik_PvdA_VVD") 
  
  coal_lik_PvdA <- gamma %>% select(coallik_CDA_PvdA, coallik_PvdA_SP_GroenLinks, coallik_PvdA_VVD) %>%
    mutate(sum_exp = exp(coallik_CDA_PvdA) + exp(coallik_PvdA_SP_GroenLinks) + exp(coallik_PvdA_VVD),
           coallik_CDA_PvdA = exp(coallik_CDA_PvdA)/sum_exp,
           coallik_PvdA_SP_GroenLinks = exp(coallik_PvdA_SP_GroenLinks)/sum_exp,
           coallik_PvdA_VVD = exp(coallik_PvdA_VVD)/sum_exp) %>%
    select(-sum_exp) %>%
    as.matrix()
  
  coal_lik_CDA <- gamma %>% select(coallik_CDA_VVD, coallik_CDA_PvdA) %>%
    mutate(sum_exp = exp(coallik_CDA_VVD) + exp(coallik_CDA_PvdA),
           coallik_CDA_VVD = exp(coallik_CDA_VVD)/sum_exp,
           coallik_CDA_PvdA = exp(coallik_CDA_PvdA)/sum_exp) %>%
    select(-sum_exp) %>%
    as.matrix()
  
  coal_lik_VVD <- gamma %>% select(coallik_CDA_VVD, coallik_PvdA_VVD) %>%
    mutate(sum_exp = exp(coallik_CDA_VVD) + exp(coallik_PvdA_VVD),
           coallik_CDA_VVD = exp(coallik_CDA_VVD)/sum_exp,
           coallik_PvdA_VVD = exp(coallik_PvdA_VVD)/sum_exp) %>%
    select(-sum_exp) %>%
    as.matrix()
  
  coal_lik_GroenLinks <- gamma %>% select(coallik_PvdA_SP_GroenLinks)  %>%
    mutate_all(as.numeric) %>%
    mutate(sum_exp =  exp(coallik_PvdA_SP_GroenLinks),
           coallik_PvdA_SP_GroenLinks = exp(coallik_PvdA_SP_GroenLinks)/sum_exp) %>%
    select(-sum_exp) %>%
    as.matrix()
  
  coal_lik_SP <- gamma %>% select(coallik_PvdA_SP_GroenLinks)  %>%
    mutate_all(as.numeric) %>%
    mutate(sum_exp =  exp(coallik_PvdA_SP_GroenLinks),
           coallik_PvdA_SP_GroenLinks = exp(coallik_PvdA_SP_GroenLinks)/sum_exp) %>%
    select(-sum_exp) %>%
    as.matrix()
  
  # Mean-variance model
  EV <- function(gamma,Z){
    gamma %*% Z
  }
  
  Vari <- function(gamma,Z){
    gamma %*% ((Z - as.numeric(EV(gamma,Z)))^2)
  }
  
  N <- nrow(dt)
  
  M_PvdA <- sapply(1:N, function(i) EV(coal_lik_PvdA[i,], rat_coal_PvdA[i,]))
  V_PvdA <- sapply(1:N, function(i) Vari(coal_lik_PvdA[i,], rat_coal_PvdA[i,]))
  
  M_CDA <- sapply(1:N, function(i) EV(coal_lik_CDA[i,], rat_coal_CDA[i,]))
  V_CDA <- sapply(1:N, function(i) Vari(coal_lik_CDA[i,], rat_coal_CDA[i,]))
  
  M_VVD <- sapply(1:N, function(i) EV(coal_lik_VVD[i,], rat_coal_VVD[i,]))
  V_VVD <- sapply(1:N, function(i) Vari(coal_lik_VVD[i,], rat_coal_VVD[i,]))
  
  M_GroenLinks <- sapply(1:N, function(i) EV(coal_lik_GroenLinks[i,], rat_coal_GroenLinks[i,]))
  V_GroenLinks <- sapply(1:N, function(i) Vari(coal_lik_GroenLinks[i,], rat_coal_GroenLinks[i,]))
  
  M_SP <- sapply(1:N, function(i) EV(coal_lik_SP[i,], rat_coal_SP[i,]))
  V_SP <- sapply(1:N, function(i) Vari(coal_lik_SP[i,], rat_coal_SP[i,]))
  
  # Gender, variablename V420
  dt$male <- ifelse(dt$V420 == "male", 1, 0)

  # Age
  dt$birthyear <- as.numeric(as.character(dt$V421)) 
  dt$age <- 2006 - dt$birthyear

  # Education
  dt$edu <- ifelse(dt$V430 == "elementary", 1, NA)
  dt$edu <- ifelse(dt$V430 == "(lower) vocational", 2, dt$edu)
  dt$edu <- ifelse(dt$V430 == "secondary", 3, dt$edu)
  dt$edu <- ifelse(dt$V430 == "middle level vocational, higher level secondary", 4, dt$edu)
  dt$edu <- ifelse(dt$V430 == "higher level vocational, university", 5, dt$edu)
  
  # PID
  if (!exists("robustnesscheck_pid")) {
    robustnesscheck_pid <- FALSE # PID as robustness?
  }
  dt$pid <- as.character(dt$V063)
  dt$pid <- ifelse(!is.na(dt$pid), dt$pid, ifelse(dt$V060=="no", "None", NA))

  dt$pid_PvdA <- ifelse(dt$pid=="PvdA", 1, 0)
  dt$pid_CDA <- ifelse(dt$pid=="CDA", 1, 0)
  dt$pid_VVD <- ifelse(dt$pid=="Partij vd Vrijheid" | dt$pid=="VVD", 1, 0)
  
  d <- data.frame(
    "vc_PvdA" = dt$vc_PvdA,
    "vc_CDA" = dt$vc_CDA,
    "vc_VVD" = dt$vc_VVD,
    "rat.PvdA" = dt$rat_PvdA,
    "rat.CDA" = dt$rat_CDA,
    "rat.VVD" = dt$rat_VVD,
    "rat.GroenLinks" = dt$rat_GroenLinks,
    "rat.SP" = dt$rat_SP,
    "pid_PvdA" = dt$pid_PvdA,
    "pid_CDA" = dt$pid_CDA,
    "pid_VVD" = dt$pid_VVD,
    "lotterymean.PvdA" = M_PvdA,
    "lotteryvariance.PvdA" = V_PvdA,
    "lotterymean.CDA" = M_CDA,
    "lotteryvariance.CDA" = V_CDA,
    "lotterymean.VVD" = M_VVD,
    "lotteryvariance.VVD" = V_VVD,
    "lotterymean.GroenLinks" = M_GroenLinks,
    "lotteryvariance.GroenLinks" = V_GroenLinks,
    "lotterymean.SP" = M_SP,
    "lotteryvariance.SP" = V_SP,
    "sex" = dt$male,
    "age" = dt$age,
    "edu" = dt$edu
  )  
  
  d$lotteryvariance.PvdA <- scales::rescale(d$lotteryvariance.PvdA, to = c(0, 1), from = c(0,0.25))
  d$lotteryvariance.CDA <- scales::rescale(d$lotteryvariance.CDA, to = c(0, 1), from = c(0,0.25))
  d$lotteryvariance.VVD <- scales::rescale(d$lotteryvariance.VVD, to = c(0, 1), from = c(0,0.25))
  
# Save model results
  summary(m1 <- lm(as.formula(paste("vc_PvdA ~ lotteryvariance.PvdA + lotterymean.PvdA + rat.PvdA + sex + as.factor(edu) + age", 
                                    if (robustnesscheck_pid) "+ pid_PvdA" else "")), d))
  summary(m2 <- lm(as.formula(paste("vc_CDA ~ lotteryvariance.CDA + lotterymean.CDA + rat.CDA + sex + as.factor(edu) + age", 
                                    if (robustnesscheck_pid) "+ pid_CDA" else "")), d))
  summary(m3 <- lm(as.formula(paste("vc_VVD ~ lotteryvariance.VVD + lotterymean.VVD + rat.VVD + sex + as.factor(edu) + age", 
                                    if (robustnesscheck_pid) "+ pid_VVD" else "")), d))
  
  NED_2006_PES_PvdA_Variance_Estimate <- tidy(m1) %>% filter(term=="lotteryvariance.PvdA") %>% select("estimate")
  NED_2006_PES_PvdA_Variance_SE <- se_robust(m1)["lotteryvariance.PvdA"]
  NED_2006_PES_PvdA_Mean_Estimate <- tidy(m1) %>% filter(term=="lotterymean.PvdA") %>% select("estimate")
  NED_2006_PES_PvdA_Mean_SE <- se_robust(m1)["lotterymean.PvdA"]
  
  NED_2006_PES_CDA_Variance_Estimate <- tidy(m2) %>% filter(term=="lotteryvariance.CDA") %>% select("estimate")
  NED_2006_PES_CDA_Variance_SE <- se_robust(m2)["lotteryvariance.CDA"]
  NED_2006_PES_CDA_Mean_Estimate <- tidy(m2) %>% filter(term=="lotterymean.CDA") %>% select("estimate")
  NED_2006_PES_CDA_Mean_SE <- se_robust(m2)["lotterymean.CDA"]
  
  NED_2006_PES_VVD_Variance_Estimate <- tidy(m3) %>% filter(term=="lotteryvariance.VVD") %>% select("estimate")
  NED_2006_PES_VVD_Variance_SE <- se_robust(m3)["lotteryvariance.VVD"]
  NED_2006_PES_VVD_Mean_Estimate <- tidy(m3) %>% filter(term=="lotterymean.VVD") %>% select("estimate")
  NED_2006_PES_VVD_Mean_SE <- se_robust(m3)["lotterymean.VVD"]
  
  # Harmonize names of the IVs of interest in the models
  names(m1$coefficients)[names(m1$coefficients) == "lotteryvariance.PvdA"] <- "Government Lottery Variance"
  names(m1$coefficients)[names(m1$coefficients) == "lotterymean.PvdA"] <- "Government Lottery Mean"
  names(m2$coefficients)[names(m2$coefficients) == "lotteryvariance.CDA"] <- "Government Lottery Variance"
  names(m2$coefficients)[names(m2$coefficients) == "lotterymean.CDA"] <- "Government Lottery Mean"
  names(m3$coefficients)[names(m3$coefficients) == "lotteryvariance.VVD"] <- "Government Lottery Variance"
  names(m3$coefficients)[names(m3$coefficients) == "lotterymean.VVD"] <- "Government Lottery Mean"
  
  
  if(!robustnesscheck_pid){
    stargazer(m1, m2, m3, title="Linear regressions of vote choice on perceived government lottery variance and mean (Dutch Parliamentary Election Study 2006)", 
              align=TRUE, 
              dep.var.labels=c("PvdA", "CDA", "VVD"),
              omit.stat=c("LL","ser","f"), 
              type = "text",
              se = lapply(list(m1, m2, m3), se_robust),
              p = lapply(list(m1, m2, m3), p_robust),
              out = "TableSM14.tex"
    )
  }
 
  